セットアップ

ライブラリーの読み込み

library(tidyverse) # recommend # see https://heavywatal.github.io/rstats/programming.html
## ─ Attaching packages ──────────────────── tidyverse 1.3.0 ─
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ─ Conflicts ───────────────────── tidyverse_conflicts() ─
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(here) # recommend # see https://uribo.hatenablog.com/entry/2018/01/25/082000
## here() starts at /Users/rakkoseminar/Dropbox/_ResearchProjects/Rproject_template_jp
library(patchwork) # recommend # see https://qiita.com/nozma/items/4512623bea296ccb74ba

関数の読み込み

# よく使う関数は、_functions_commonフォルダに作成
# プロジェクトで使う関数は、_functions_projectフォルダに作成し、source()で読み込むこと
common_functions_path <- here("Analysis_Script", "_functions_common")
common_functions <- here(common_functions_path, dir(common_functions_path))
purrr::walk(common_functions, ~source(.x))

project_functions_path <- here("Analysis_Script", "_functions_project")
project_functions <- here(project_functions_path, dir(project_functions_path))
purrr::walk(project_functions, ~source(.x))
# source()は複数の関数を読み込めないので、purrr::walk()を使う
# purrr::walk()は、vector (or list) を関数に適用しながらinvisible返しする、purrr::mapの亜種

メタ情報の設定

# knit時の各種設定
# 特に、cacheを _cacheフォルダに生成するための設定
# このchunkは手動でrunしないこと。エラーが出る。

rmarkdown_dir <- "Analysis_Script"
cache_dir <- "_cache" # cacheの保存ディレクトリ名
input <- knitr::current_input() %>% str_split(pattern = "\\.Rmd") %>% .[[1]] %>% .[1] # 拡張子の前だけ取り出す # 返り値がリストになっているのでその中身を取り出す

knitr::opts_chunk$set(
  echo=TRUE, comment="", warning=FALSE, message=FALSE, fig.align="center",
  
  cache.path = here(rmarkdown_dir, cache_dir, input %&% "_cache", "html") %&% "/", 
  fig.path = here(rmarkdown_dir, cache_dir, input %&% "_files", "figure-html") %&% "/"
  )
# 出力の保存先のパスを設定
# dir.create(here("Output", "Figures_DummyData")) # なければ、フォルダを作成
output_path <- here("Output", "Figures_DummyData")

データの読み込み

データ読み込みの基本的な方針

  • このRmarkdown内で使用するデータは、なるべくこのセクション内で全て読み込み、データフレームとして保存しておくこと

  • もしRawDataがあまり整理されていない場合には、整形用のRスクリプトを別途用意して、Rの中でデータを整形すること。そして、その結果を新しいcsvファイルとして ./Data/に保存し、それを読み込み直すのが良い。整形コードは、Rmarkdownに書かない方が良い。

  • 2種類以上のデータ分析を行い、それぞれ全く別のデータフレームが分析対象となるとしても、このセクションでまとめて読み込み、別々のデータフレームとして保存しておくのが良い。

    • あまりおすすめしない方法:分析1の直前に分析1で使用するデータを読み込み、分析2の直前に分析2で使用するデータを読みこむ。
    • 理由:ここら辺は好みにもよるが、あとで見返した時にデータフレームの読み込みがどこで行われたがわかりにくくなりがち。
    • ただし、あまりにもデータ量が多いため、複数のオブジェクトを保持したくない、という場合もある(毎回同じオブジェクト名を使用して、中身だけ更新したい)。これはケースバイケースで対応する。

※ 分析の際には、上記の文章を削除すること

# --- 読み取りの例 --- # 
# デフォルトで、Data/DummyData_RL/に3つのダミーデータのディレクトリがあるはず。  
# テストとして、それらを読み込んでデータフレームにする。  
# ディレクトリがなければ、generate_DummyData_RL.R or generate_DummyData_RL.py のどちらかを実行してcsvを生成しておく。

options(readr.num_columns = 0) # readrのメッセージがうるさいので消す

DataDir <- "Data"
DataDir_sub <- "DummyData_RL"

target_csv <- make_DataPath(DataDir, DataDir_sub, ".", ".csv")
df_dummy_RL <- purrr::map_dfr(target_csv, readr::read_csv)
# 複数のcsv(列構造が同じもの)は、purrr::map_dfr()がrapid & beautifulで良い
# csvが巨大な場合には、purrr::map_dfr(target_csv, data.table::fread)を実行すると良い
options(readr.num_columns = 0) # readrのメッセージがうるさいので消す

df_est_params <- readr::read_csv(here("Data", "Models_Results", "basic_RL", "fit01", "estimated_parameter", "est_parameter.csv")) %>% 
  dplyr::rename(
    ID = p_index
  ) %>% 
  dplyr::filter(para %in% c("alpha", "beta")) %>% 
  tidyr::pivot_wider(id_cols = c("ID"), names_from = para, values_from = mean) %>% 
  dplyr::rename(
    fit_alpha = alpha,
    fit_beta = beta
  )

概要

これまでの分析の要約、今回の分析の目的、どんな分析をやったのかの要約を列挙する

分析1:試行を通してのQ値と選択の変化

画像の作成

作図の基本的な方針

  • 横軸や縦軸などのパラメータが同じで、与えるデータフレームだけが異なる場合は、作図を関数化して、purrr::map()でまとめて作成する。または、 %+% でデータだけ置き換えることで作成する。
  • 論文や学会スライドの際にすぐ利用できるように、なるべくggsave()でpng形式で保存しておく。ファイル名には一定の規則を設けて出力する。
  • Figureの保存先は、メタ情報の設定:setup-output-pathチャンクで指定済み。

※ 分析の際には、上記の文章を削除すること

# 横軸:試行、縦軸:Q値
## 1人分の画像なので、purrr::map()で各課題、各参加者分の画像をまとめて生成すること

library(tidyverse)

plot_TrialQvalue <- function(df) {
  
  g <- df %>% 
    ggplot(aes(x = trial, y = value, color = as.factor(Q))) + 
    geom_line(size = 1.5) + 
    my_theme2 + 
    theme(
      axis.text.x = element_text(size = 18), 
      axis.text.y = element_text(size = 15), 
      strip.text.x = element_text(size = 15),
      strip.text.y = element_text(size = 15)) +
    scale_x_continuous(breaks = c(0, unique(df$trial)[length(unique(df$trial))]/2, unique(df$trial)[length(unique(df$trial))])) + 
    coord_cartesian(ylim = c(0, 1.0)) + 
    labs(x = "trial", y = "Q value", color = "Q") 
  
  return(g)
}
# 横軸:試行、縦軸:Q値

#dir.create(here(output_path, "Qvalue"))

## 課題ごと(報酬確率が異なる)、facetで行:beta, 列:alphaにする

task_sort <- c("Op0.6", "Op0.7", "Op0.8")

df_dummy_RL %>% 
  tidyr::pivot_longer(cols = c(Q1, Q2), names_to = "Q", values_to = "value") %>% 
  dplyr::mutate(
    beta_c = "beta = " %&% beta,
    alpha_c = "alpha = " %&% alpha,
  ) %>% 
  dplyr::group_split(Op1_p) -> 
  df_dummy_RL_list 

# 画像の作成
g <- purrr::map(df_dummy_RL_list, 
                ~{p <- plot_TrialQvalue(.x); 
                  p + facet_grid(beta_c ~ alpha_c)})

# 画像の保存
purrr::walk2(task_sort, g, 
             ~ggsave(filename = here(output_path, "Qvalue", "Qvalue_full_" %&% .x %&% ".png"),
                     plot = .y, height = 16, width = 18))

## 課題ごと(報酬確率が異なる)、IDごと(パラメータが異なる)に分割してlistにし、purrr::map()でまとめて画像を生成

df_dummy_RL %>% 
  tidyr::pivot_longer(cols = c(Q1, Q2), names_to = "Q", values_to = "value") %>% 
  dplyr::group_split(Op1_p, ID) ->
  df_dummy_RL_list

IDs <- 1:25
plot_params <- # 保存用のdf
  tibble(
    task = task_sort
  ) %>% 
  tidyr::crossing(ID = IDs) %>% 
  dplyr::mutate(
    alpha = df_dummy_RL_list %>% purrr::map("alpha") %>% purrr::map_dbl(1),
    beta = df_dummy_RL_list %>% purrr::map("beta") %>% purrr::map_dbl(1)
  )
  
# 画像の作成
g <- purrr::map(1:nrow(plot_params), 
                ~plot_TrialQvalue(df = df_dummy_RL_list[[.x]]))

# 画像の保存
purrr::walk2(1:nrow(plot_params), g, 
             ~{save_name <- "Qvalue_task" %&% plot_params[.x,]$task %&% "_ID" %&% plot_params[.x,]$ID %&% "_alpha" %&% plot_params[.x,]$alpha %&% "_beta" %&% plot_params[.x,]$beta %&% ".png"; 
             ggsave(filename = here(output_path, "Qvalue", save_name),
                     plot = .y, height = 6, width = 8)})
# 横軸:試行、縦軸:選択確率と選択
## 1人分の画像なので、purrr::map()で各課題、各参加者分の画像をまとめて生成すること

library(tidyverse)

plot_TrialProb <- function(df) {
  
  g <- df %>% 
    dplyr::mutate(
      Choice = if_else(Choice == 2, 0, 1) # choice 2 -> 0
    ) %>% 
    ggplot() + 
    geom_point(aes(x = trial, y = Choice, shape = as.factor(Reward))) + 
    geom_line(aes(x = trial, y = Prob), size = 1.5) + 
    my_theme2 + 
    theme(
      axis.text.x = element_text(size = 18), 
      axis.text.y = element_text(size = 15), 
      strip.text.x = element_text(size = 15),
      strip.text.y = element_text(size = 15)) +
    scale_x_continuous(breaks = c(0, unique(df$trial)[length(unique(df$trial))]/2, unique(df$trial)[length(unique(df$trial))])) + 
    coord_cartesian(ylim = c(0, 1.0)) + 
    labs(x = "trial", y = "Probability of choosing 1", shape = "Reward") 
  
  return(g)
  
}
# 横軸:試行、縦軸:p値と選択

dir.create(here(output_path, "CProb"))

## 課題ごと(報酬確率が異なる)、facetで行:beta, 列:alphaにする

task_sort <- c("Op0.6", "Op0.7", "Op0.8")

df_dummy_RL %>% 
  dplyr::mutate(
    beta_c = "beta = " %&% beta,
    alpha_c = "alpha = " %&% alpha,
  ) %>% 
  dplyr::group_split(Op1_p) -> 
  df_dummy_RL_list 

# 画像の作成
g <- purrr::map(df_dummy_RL_list, 
                ~{p <- plot_TrialProb(.x); 
                  p + facet_grid(beta_c ~ alpha_c)})

# 画像の保存
purrr::walk2(task_sort, g, 
             ~ggsave(filename = here(output_path, "CProb", "CProb_full_" %&% .x %&% ".png"),
                     plot = .y, height = 16, width = 18))

## 課題ごと(報酬確率が異なる)、IDごと(パラメータが異なる)に分割してlistにし、purrr::map()でまとめて画像を生成

df_dummy_RL %>% 
  dplyr::group_split(Op1_p, ID) ->
  df_dummy_RL_list

IDs <- 1:25
plot_params <- # 保存用のdf
  tibble(
    task = task_sort
  ) %>% 
  tidyr::crossing(ID = IDs) %>% 
  dplyr::mutate(
    alpha = df_dummy_RL_list %>% purrr::map("alpha") %>% purrr::map_dbl(1),
    beta = df_dummy_RL_list %>% purrr::map("beta") %>% purrr::map_dbl(1)
  )
  
# 画像の作成
g <- purrr::map(1:nrow(plot_params), 
                ~plot_TrialProb(df = df_dummy_RL_list[[.x]]))

# 画像の保存
purrr::walk2(1:nrow(plot_params), g, 
             ~{save_name <- "CProb_task" %&% plot_params[.x,]$task %&% "_ID" %&% plot_params[.x,]$ID %&% "_alpha" %&% plot_params[.x,]$alpha %&% "_beta" %&% plot_params[.x,]$beta %&% ".png"; 
             ggsave(filename = here(output_path, "CProb", save_name),
                     plot = .y, height = 6, width = 8)})

報酬確率0.8 vs 0.2

報酬確率0.7 vs 0.3

報酬確率0.6 vs 0.4

分析2:パラメータの推定結果と真値の相関

  • 横軸:シミュレーション時のパラメータの値
  • 縦軸:推定したパラメータの値(平均値)
# データの処理
df_dummy_RL %>% 
  dplyr::filter(Op1_p == 0.8) %>% 
  dplyr::filter(trial == 1) %>% 
  dplyr::transmute(
    ID = ID,
    sim_alpha = alpha,
    sim_beta = beta
  ) %>% 
  dplyr::full_join(df_est_params, by = "ID") -> 
  df_sim_fit
df_sim_fit %>% 
  ggplot(aes(x = sim_alpha, y = fit_alpha)) + 
  geom_point() + geom_abline(slope = 1, intercept = 0) + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1))

df_sim_fit %>% 
  ggplot(aes(x = sim_beta, y = fit_beta)) + 
  geom_point() + geom_abline(slope = 1, intercept = 0)

SessionInfo

sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] patchwork_1.1.1 here_1.0.1      forcats_0.5.0   stringr_1.4.0  
 [5] dplyr_1.0.2     purrr_0.3.4     readr_1.4.0     tidyr_1.1.2    
 [9] tibble_3.0.4    ggplot2_3.3.2   tidyverse_1.3.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.0  xfun_0.19         haven_2.3.1       colorspace_2.0-0 
 [5] vctrs_0.3.5       generics_0.1.0    htmltools_0.5.1.1 yaml_2.2.1       
 [9] rlang_0.4.9       pillar_1.4.7      glue_1.4.2        withr_2.3.0      
[13] DBI_1.1.0         dbplyr_2.0.0      modelr_0.1.8      readxl_1.3.1     
[17] lifecycle_0.2.0   munsell_0.5.0     gtable_0.3.0      cellranger_1.1.0 
[21] rvest_0.3.6       evaluate_0.14     labeling_0.4.2    knitr_1.30       
[25] fansi_0.4.1       broom_0.7.2       Rcpp_1.0.5        scales_1.1.1     
[29] backports_1.2.0   jsonlite_1.7.1    farver_2.0.3      fs_1.5.0         
[33] hms_0.5.3         digest_0.6.27     stringi_1.5.3     grid_4.0.2       
[37] rprojroot_2.0.2   cli_2.2.0         tools_4.0.2       magrittr_2.0.1   
[41] crayon_1.3.4      pkgconfig_2.0.3   ellipsis_0.3.1    xml2_1.3.2       
[45] reprex_0.3.0      lubridate_1.7.9.2 assertthat_0.2.1  rmarkdown_2.5    
[49] httr_1.4.2        rstudioapi_0.13   R6_2.5.0          compiler_4.0.2